1 Wstęp

Surowy zbiór danych był na tyle duży, że problematyczne było załadowanie go w całości do pamięci RAM. Dopiero usunięcie pewnej części zbędnych danych jeszcze przed wczytaniem pozwoliło na wczytanie go do data frame’a. Mimo usunięcia sporej liczby krotek i atrybutów podczas wstępnego przetwarzania, nadal w fazie uczenia maszynowego mieliśmy doczynienia z ponad 350000 obserwacjami, po 383 atrybuty każda. Było to zbyt dużo dla większości algorytmów, aby być w stanie zbudować model korzystając z pozostałej wolnej pamięci RAM. W tym wypadku pierwotnie zdecydowano się na podejście przyrostowe i budowanie modelu stopniowo z na tyle małych wycinków zbioru danych, aby było możliwe przeprowadzenie uczenia. Później jednak wybrano selekcję atrybutów jak sposób na rozwiązanie problemów z pamięcią. Przy okazji przyspieszyło to znacznie same obliczenia.

2 Ładowanie bibliotek

library(dplyr)
library(ggplot2)
library(caret)
library(DT)
library(reshape2)
library(tidyr)
library(plotly)
library(biglm)
library(xgboost)

3 Wstępne przetwarzanie

set.seed(95)

3.1 Wczytywanie i czyszczenie danych

Usunięcie atrybutów redundantnych, pustych, o zbyt złożonej strukturze do analizy. Decyzję podjęto na podstawie wstępnego przejrzenia pliku all_summary.csv zawierającego cały zbiór danych. Pozwala to zmniejszyć rozmiar wczytywanego pliku, a w konsekwencji rozmiar wynikowego data frame’a. Dzięki temu zabiegowi jest on w stanie zmieścić się w pamięci RAM.

if (!("data.csv" %in% dir())) 
    system("cut -d';' -f1-3,27,399 --complement all_summary/all_summary.csv > data.csv")

Wczytanie danych i usunięcie kolumn które nie wezmą udziału w dalszym przetwarzaniu.

df <- read.table("data.csv", nrows = 75000, header = T, sep=";", comment.char = "")
classes <- sapply(df, class)
df <- read.table("data.csv", nrows = 591043, header = T, sep=";", comment.char = "", colClasses = classes)
colToRemove <- c("local_BAa", "local_NPa", "local_Ra", "local_RGa", "local_SRGa", "local_CCSa", "local_CCPa", "local_ZOa", "local_ZDa", "local_ZD_minus_a", "local_ZD_plus_a", "local_res_atom_count", "local_res_atom_non_h_occupancy_sum", "local_res_atom_non_h_electron_occupancy_sum", "local_res_atom_C_count", "local_res_atom_N_count", "local_res_atom_O_count", "local_res_atom_S_count", "dict_atom_C_count", "dict_atom_N_count", "dict_atom_O_count", "dict_atom_S_count", "fo_col", "fc_col", "weight_col", "grid_space", "solvent_radius", "solvent_opening_radius", "part_step_FoFc_std_min", "part_step_FoFc_std_max", "part_step_FoFc_std_step")
df <- df[, !(names(df) %in% colToRemove)]

Usunięcie ze zbioru danych ktrotek, które posiadają wskazane wartości atrybutu res_name

df <- filter(df, !res_name %in% c("UNK", "UNX", "UNL", "DUM", "N", "BLOB", "ALA", "ARG", "ASN", "ASP", "CYS", "GLN", "GLU", "GLY", "HIS", "ILE", "LEU", "LYS", "MET", "MSE", "PHE", "PRO", "SEC", "SER", "THR", "TRP", "TYR", "VAL", "DA", "DG", "DT", "DC", "DU", "A", "G", "T", "C", "U", "HOH", "H20", "WAT"))

3.2 Obsługa wartości pustych

Na tym etapie przetwarzania w zbiorze danych występuje 40036 krotrek o conajmniej jednej wartości pustej. Przekłada się to na nieco mniej niż 7% całego zbioru. Ponieważ dalej chcemy rozważać tylko 50 najczęstszych wartości atrybutu res_name sprawdzimy, czy usunięcie wszytskich krotek z wartościami pustymi wpłynie na ten wybór. Okazuje się, że po zostawieniu jedynie kompletnych krotek nadal wybierzemy te same 50 wartości res_name. W związku z tym decydujemy się na tę obsługę wartości pustych.

sumna <- sapply(df, function(x) sum(is.na(x)))
sumna <- sumna[sumna > 0]
print(max(sumna))
## [1] 40036
print(max(sumna)/length(df[[1]]))
## [1] 0.06839797
dfComplete <- df[complete.cases(df),]
countn <- (df %>% group_by(res_name) %>% summarise(n()) %>% rename(sum = "n()") %>% filter(!is.na(res_name)) %>% arrange(desc(sum)))[1:50,1]
countn2 <- (dfComplete %>% group_by(res_name) %>% summarise(n()) %>% rename(sum = "n()") %>% arrange(desc(sum)))[1:50,1]
length(intersect(countn$res_name, countn2$res_name))
## [1] 50
df <- dfComplete
rm(dfComplete)

3.3 Podsumowanie przetworzonych danych

Po wstępnym przetworzeniu zbiór danych składa się z 525666 krotek, każda zawiera 388 atrybutów. 4 z nich są typu factor a pozostałe 384 to liczby, gdzie 51 to liczby całkowite a 333 zmiennoprzecinkowe.

print(length(df[[1]]))
## [1] 525666
print(length(df))
## [1] 388
knitr::kable(table(sapply(df, class)))
Var1 Freq
factor 4
integer 51
numeric 333

4 Analiza danych

Sprawdzenie liczności każdej z 50 klas.

df <- filter(df, res_name %in% countn2$res_name)
countn <- df %>% group_by(res_name) %>% summarise(n()) %>% rename(sum = "n()") %>% arrange(desc(sum))
df$res_name <- factor(df$res_name)
prettyTable(countn)

4.1 Korelacja

Po wyliczeniu korelacji między atrybutami liczbowymi okazuje się, że wiele z nich jest bardzo silnie skorelowanych.

correl <- cor(df[, sapply(df, is.numeric)])
correl[lower.tri(correl)] <- NA
correl <- melt(correl)
correl <- correl[complete.cases(correl), ] %>% filter(Var1 != Var2) %>% mutate(absValue = abs(value)) %>% arrange(desc(absValue))
prettyTable(correl)
cantUse <- c("dict_atom_non_h_count", "dict_atom_non_h_electron_sum", "local_res_atom_non_h_electron_sum", "local_res_atom_non_h_count")
colForElectron <- (correl %>% filter(Var1 == "local_res_atom_non_h_count", abs(value) > 0.6))$Var2
colForElectron <- as.vector(colForElectron[!colForElectron %in% cantUse])
colForAtom <- (correl %>% filter(Var1 == "local_res_atom_non_h_electron_sum", abs(value) > 0.6))$Var2
colForAtom <- as.vector(colForAtom[!colForAtom %in% cantUse])
rm(correl)

4.2 Wykresy rozkładów

Zarówno rozkład liczby atomów jak i elektronów skupia się głównie na mniejszych wartościach, stosunkowo niewielkich w porównaniu do zakresu wartości na tych atrybutach.

ggplotly(ggplot(df, aes(local_res_atom_non_h_count)) + geom_density())
ggplotly(ggplot(df, aes(local_res_atom_non_h_electron_sum)) + geom_density())

4.3 Analiza niezgodności liczby atomów i elektronów

Ze względu na stosunkowo duży zakres wartości obu atrybutów, zdecydowano się bazować na błędzie względnym zamiast bezwzględnego. Podejście te ma tą własność, że pomyłka o 2 elektrony będzie dużym błędem gdy rzeczywista wartość to 1, natomiast znacznie mniejszym przy rzeczywistej wartości 400.

knitr::kable((df %>% group_by(res_name) %>% summarise(mean(abs(1-local_res_atom_non_h_count/dict_atom_non_h_count))) %>% rename(sum = 2) %>% arrange(desc(sum)))[1:10,] %>% rename("średni błąd względny liczby atomów" = 2))
res_name średni błąd względny liczby atomów
1PE 0.1671000
MLY 0.1093553
SEP 0.0909731
PG4 0.0799353
MAN 0.0723490
NAG 0.0653264
CLA 0.0550504
PLP 0.0539990
PGE 0.0402994
COA 0.0392920
knitr::kable((df %>% group_by(res_name) %>% summarise(mean(abs(1-local_res_atom_non_h_electron_sum/dict_atom_non_h_electron_sum))) %>% rename(sum = 2) %>% arrange(desc(sum)))[1:10,] %>% rename("średni błąd względny liczby elektronów" = 2))
res_name średni błąd względny liczby elektronów
1PE 0.1688296
MLY 0.1267699
SEP 0.0914773
MAN 0.0826846
PG4 0.0798827
NAG 0.0760858
PLP 0.0583479
CLA 0.0527664
PGE 0.0396442
COA 0.0389528

4.3 Rozkałd wartości atrybutów part_01

5 Uczenie maszynowe

W tej sekcji nastąpi próba obliczania liczby atomów i elektronów, oraz przewidzenia wartości atrybutu res_id. ##Regresja liczby atomów

moreColToRemove <- c("dict_atom_non_h_count", "dict_atom_non_h_electron_sum", "pdb_code", "res_id", "chain_id")
df <- df[, !(names(df) %in% moreColToRemove)]
forTrain <- createDataPartition(df$local_res_atom_non_h_count, p= .95, list=F)
formula <- as.formula(paste("local_res_atom_non_h_count ~ ", paste(colForAtom, collapse = " + ")))
split = 20
colForAtom <- append(colForAtom, "local_res_atom_non_h_count")
df2 <- df[, names(df) %in% colForAtom]
model <- lm(formula, data=df2[forTrain, ])
y <- predict(model, df2[-forTrain, ])
## Warning in predict.lm(model, df2[-forTrain, ]): prediction from a rank-
## deficient fit may be misleading
diff <- y-df[-forTrain, "local_res_atom_non_h_count"]
diff <- diff*diff
rmse <- sqrt(mean(diff))
r2 <- cor(y, df[-forTrain, "local_res_atom_non_h_count"])^2
print(rmse)
## [1] 9.382383
print(r2)
## [1] 0.4902731
rm(model, y, diff)

5.1 Regresja liczby elektronów

forTrain <- createDataPartition(df$local_res_atom_non_h_electron_sum, p= .95, list=F)
formula <- as.formula(paste("local_res_atom_non_h_electron_sum ~ ", paste(colForElectron, collapse = " + ")))
split = 20
colForElectron <- append(colForElectron, "local_res_atom_non_h_electron_sum")
df2 <- df[, names(df) %in% colForElectron]
model <- biglm(formula, df2[forTrain, ])
y <- predict(model, df2[-forTrain, ])
diff <- y-df[-forTrain, "local_res_atom_non_h_electron_sum"]
diff <- diff*diff
rmse <- sqrt(mean(diff))
r2 <- cor(y, df[-forTrain, "local_res_atom_non_h_electron_sum"])^2
print(rmse)
## [1] 65.45624
print(r2)
##           [,1]
## [1,] 0.4620111
rm(model, y, diff)

Regresja nie okazała się być bardzo dokładną, jednak może to być spowodowane prostotą modelu linowego. ##Klasyfikacja W celu poprawienia jakości klasyfikacji, a przede wszystkim czasu uczenia i zajętości pamięciowej dokonana została selekcja atrybutów. Najpierw usunięto atrybuty o zerowej wariancji a następnie te, które były silnie skorelowane z innymi. Jako klasyfikatora użyto xgboost.

colToRemove <- c("local_res_atom_non_h_count", "local_res_atom_non_h_electron_sum")
df <- df[, !(names(df) %in% colToRemove)]
df <- df[,sapply(df, function(v) !is.numeric(v) || (var(v)!=0))]
correl <- cor(df[, sapply(df, is.numeric)])
correl[lower.tri(correl, diag = T)] <- 0
correl <- apply(correl, 2, max)
correl <- correl[correl > 0.8]
df <- df[, !(names(df) %in% names(correl))]
df$res_name <- factor(df$res_name)
forTrain <- createDataPartition(df$res_name, p=.7, list = F)
dtest <- xgb.DMatrix(data = as.matrix(df[-forTrain, names(df) != "res_name"]), label = df[-forTrain, "res_name"])
dtrain <- xgb.DMatrix(data = as.matrix(df[forTrain, names(df) != "res_name"]), label = df[forTrain, "res_name"])
watchlist <- list(train=dtrain, test=dtest)
bstDMatrix <- xgb.train(data = dtrain, max.depth = 10, eta = 1, nthread = 4, nrounds = 1, objective = "multi:softmax", tree_method="approx", num_class=51, watchlist=watchlist)
## [23:50:39] Tree method is selected to be 'approx'
## [1]  train-merror:0.359822   test-merror:0.400867